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We compare a newly developed hybrid simulation method which combines classical molecular dy- 
namics (MD) and computational fluid dynamics (CFD) to a simulation consisting only of molecular 
dynamics. The hybrid code is composed of three regions: a classical MD region, a continuum domain 
where the dynamical equations are solved by standard CFD methods, and an overlap domain where 
transport information from the other two domains is exchanged. The exchange of information in the 
overlap region ensures that momentum, energy and mass are conserved. The validity of the hybrid 
code is demonstrated by studying a single polymer tethered to a hard wall immersed in explicit 
solvent and undergoing shear flow. In classical molecular dynamics simulation a great deal of com- 
putational time is devoted to simulating solvent molecules, although the solvent itself is of no direct 
interest. By contrast, the hybrid code simulates the polymer and surrounding solvent explicitly, 
whereas the solvent farther away from the polymer is modeled using a continuum description. In 
the hybrid simulations the MD domain is an open system whose number of particles is controlled 
to filter the perturbative density waves produced by the polymer motion. We compare conforma- 
tional properties of the polymer in both simulations for various shear rates. In all cases polymer 
properties compare extremely well between the two simulation scenarios, thereby demonstrating 
that this hybrid method is a useful way to model a system with polymers and under nonzero flow 
conditions. There is also good agreement between the MD and hybrid schemes and experimental 
data on tethered DNA in flow. The computational cost of the hybrid protocol can be reduced to 
less than 6% of the cost of updating the MD forces, confirming the practical value of the method. 

PACS numbers: 02.07.Ns, 68.05. Cf 



I. INTRODUCTION 

Molecular dynamics (MD) simulations have long been 
used to model complex fluids both in and out of equilib- 
rium. As computers get more powerful there has been an 
increasing desire for more chemically accurate models of 
these fluids. This means that simulations are becoming 
larger and more accurate, but also that much simula- 
tion time is being devoted to model in detail parts of 
the computational system of little direct scientific inter- 
est. Hybrid methods combine regions of relatively high 
degree of chemical accuracy in a specific region of in- 
terest, and a more coarse-grained model where the dy- 
namics can be solved in a less computationally intensive 
way, further away from the specific region of interest. 
We focus, in particular, on hybrid methods that combine 
Lennard-Jones type specificity with larger scale contin- 
uum methods. Such hybrid methods have been applied 
in a number of fields, including Lennard-Jones fluids l^l, 
biophysics 2] and MD/CFD couphng 3, 4J. This type 
of simulation technique is particularly useful in study- 
ing interface problems, where the region of interest is a 
localized part of the entire system. 

Typical hybrid methods consist of three regions: a tra- 
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ditional region where dynamics are simulated using well- 
established techniques such as molecular dynamics 0, 
a continuum region where CFD or elasticity differential 
equations are solved using classical techniques, and an 
overlap region where the necessary transport information 
of the MD and continuum regions are exchanged. The 
primary motivation for using a hybrid scheme is to reduce 
computer time devoted to simulating bulk regions of lit- 
tle direct interest. As such, a hybrid scheme is ideally 
suited to studying interfacial systems. 

In this paper we apply a hybrid technique to a sin- 
gle polymer tethered to a wall with explicit solvent. 
The complex dynamics arising from this system have at- 
tracted a degree of interest from experimentalists, who 
used fluorescence microscopy and videomicroscopy to 
investigate the dynamic properties of individual DNA 
chains in a shear flow, either tethered to a wall jfl] or 
free 01 ■ These experiments reveal that the structural 
quantities, such as the mean elongation of the polymer, 
are very sensitive to flow environment and that the dy- 
namical properties depend strongly on the initial con- 
formation. Moreover, care needs to be taken to control 
the finite size effects, such as those due to long-ranged 
hydrodynamic interaction between the polymer and the 
walls |7J . This large "sensitivity" of the tethered polymer 
dynamics is in fact a valuable test for the hybrid model. 
First, the hybrid model reduces the size of the MD sim- 
ulation box while avoiding finite size effects and, second, 
the coupling has to be able to perfectly reproduce flows 
at very small shear rates. As shown in recent work |^, 
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this second task is non-trivial because the signal-to-noise 
ratio of the stress that one needs to communicate from 
the particle to the continuum system is very small. 

The problem of tethered polymers under flow has a 
geometry which is ideal for a hybrid scheme. The sci- 
entific interest lies around the polymer although, in a 
standard MD simulation, the solvent particles within the 
bulk flow require most of the computational time. Sin- 
gle polymers in a bath of explicit molecular solvent have 
been the focus of a grea t deal of attention in the last 
decade 0, IH O 111 d . Many of these studies are 
devoted to examining a free chain in solution in order to 
make comparisons with theoretical predictions, explore 
the dynamics regime beyond the short-time dynamics or 
extract scaling laws as a function of polymer length or 
shear flow. In these studies, the solvent is explicity sim- 
ulated. For example, the study by Diinweg and Kremer 
^0] uses a polymer of length L = 60 beads in a bath 
of 7940 Lennard-Jones spheres. Aust, Kroger and Hess 
[m simulate polymers of length L = 10 to 60 in sys- 
tems where the total number of particles including sol- 
vent ranges from 1000 to 5832. It is clear in these cases 
that most of the computational effort is devoted to solv- 
ing the equations of motion of the solvent particles when 
the real scientific interest lies in the polymer behavior. 

The single polymer we study is tethered to a wall, and 
a variety of shear rates is imposed as a model interfacial 
problem to compare classical MD techniques to the hy- 
brid simulation. In classical MD, we sandwich the poly- 
mer and solution between two explicit walls, and impose 
periodic boundary conditions in the remaining two di- 
rections. The polymer is tethered to the bottom wall, 
and shear is created by moving the top wall at constant 
velocity in a direction parallel to the wall. In the hy- 
brid case, we model one wall, the polymer and some of 
the solvent explicitly using MD, and impose shear by a 
boundary condition in the CFD regime of the calculation. 
The shear is translated to the MD regime via energy and 
momentum flux transfers in the overlap region. We com- 
pare various conformational properties of the polymer for 
the two techniques. 

Our paper is organized as follows. In the following sec- 
tion we briefly outline both the classical MD simulation 
and the hybrid simulation techniques. In Sec. IIIII we 
compare the conformation of the polymer as calculated 
by each simulation method. The computational costs and 
benefits of the hybrid scheme are compared to classical 
MD. We also compare our results to experimental data 
of tethered DNA under shear flow. We conclude with a 
discussion in Sec. IIVI 



METHOD 



II. 



We describe in this section both the molecular dynam- 
ics and hybrid dynamics models used in our simulations. 
The MD part of the hybrid scheme was the same as the 
classical MD used in the pure molecular dynamics simu- 



lations. 



1. Molecular Dynamics 

The polymer model and simulation techniques are sim- 
ilar to those used in previous work [Ulil. The polymer 
potential is based on the bead-spring model developed 
by Kremer and Grest |17| |. Linear polymers containing 

= 60 beads each are created by linking nearest neigh- 
bors on a chain with the potential 
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where is the distance between beads i and j, Ro = 
1.5(7, k = 30e/CT^, and a and e set the length and energy 
scales, respectively. The monomers in the solvent and in 
the polymer interact through a truncated Lennard-Jones 
(LJ) potential 



(2) 



The cutoff is set at Tc = for all fluid particles to 

produce a purely repulsive interaction between beads. 

The bounds of the simulation cell are periodic in the 
X and y directions, with periods = 38.5cr and Ly ^ 
33.4(7, respectively. In the z direction the cell is bounded 
by top and bottom walls. Each wall contains two layers 
of 1600 spheres strongly tethered to the sites of a (1, 1, 1) 
plane of an fee lattice by harmonic springs of stiffness 
K = 1320e(7^^. The wall atoms do not interact with 
each other, and the wall- fluid interaction is LJ with an 
increased cutoff of Tc = 1.25(7 and increased energy scale 
of e^,/ = vTje. The increased cutoff and energy ensure 
sufficient adhesion of the fluid to the wall so that the slip 
at the wall is minimized for the shear rates considered 
here. The polymer is anchored to the wall by enforcing 
the tethering potential, Eq. between the end of the 
polymer and one wall atom. 

The walls are 48(7 apart for the pure MD simulation; 
in the hybrid simulation the molecular dynamics region 
persists for 19(7. There are sufficient solvent monomers to 
yield a mean fluid density of approximately p = 0.8(7~^ 
in the center of the simulation cell, although density os- 
cillations are induced within a few a of the walls [13 ■ 

The equations of motion are integrated using a velocity 
Verlet algorithm 5], with a time step 5t = 0.0075t, where 
T = a ^fmji. is the basic unit of time, and m is the mass 
of a monomer. A constant temperature of /cbT = l.Oe is 
maintained with a Langevin thermostat |17|. To ensure 
that this thermostat does not bias the shear profile, the 
Gaussian white noise and damping terms are only added 
to the equations of motion for the velocity components 
normal to the mean flow, that is the y and z directions 

m. 
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The shear flow in our pure MD simulation is induced by 
moving the atomistic top wall at a constant speed in 
the X direction. In the hybrid simulation, a shear bound- 
ary condition is used for the continuum regime, and this 
resulted in shear in the MD regime by the exchange of 
momentum in the overlap domain. The starting configu- 
ration was that of a single polymer tethered to the wall, 
in an equilibriated solvent. We repeat the simulation for 
two different starting configurations, i.e. each configura- 
tion has a polymer tethered to a different wall atom, and 
the initial conformation of the polymer is different. The 
initial polymer configurations were either taken from pre- 
vious simulations on melts |l6j | , or generated from a ran- 
dom walk. Although over long periods of time we expect 
that different starting configurations will give the same 
configurational averages, previous work jlQj has shown 
that hundreds of different initial configurations are re- 
quired to arrive at reasonable ensemble averages. In view 
of this, we used two initial configurations for the pure 
MD simulations; although this falls short of the number 
of initial configurations required to achieve ensemble av- 
erages, it does give us a window over which to compare 
the hybrid simulation. 

The local shear rate, 7, of the fluid is calculated by 
computing the local change in the x component of veloc- 
ity, as a function of z, i.e. 7 = dvx/dz. The upper 
wall velocity was chosen so that the shear rate 7 assumed 
the values 0.0, 0.0005, 0.001, 0.002, 0.005, O.Olr-i. Sim- 
ulations at higher shear rates were created from lower 
shear rates by increasing the wall velocity or bound- 
ary condition and allowing the system to achieve steady 
state. The simulations were done for at least one million 
time steps, and the runs of 7 = O.OOIt"^ and 0.005r^^ 
were simulated for at least ten million time steps, corre- 
sponding to a total run time of 75, OOOr. In the analysis 
in the following section the first 250, 000 time steps of 
data for each given shear rate were ignored, to allow the 
system to reach steady state; this length of time was de- 
termined to be the longest time necessary for the system 
to reach steady state once a new shear rate was imposed. 



2. Hybrid Dynamics 

Our hybrid dynamics code consists of three do- 
mains: the particle domain (P) which was by the same 
molecular dynamics method described above, the con- 
tinuum domain (C) treated by standard continuum fluid 
dynamics and, an overlap region where information from 
the other two domains is exchanged. The hybrid scheme 
is a protocol to exchange fluxes of conserved quantities, 
specifically mass, momentum and energy between both 
classically treated regimes. To implement the two-way 
flux exchange, the overlap region consists of two differ- 
ent subdomains: the P^C and the C^P cells. Within 
the C— >P region, the fluxes from the continuum domain 
are imposed on the particle domain, whereas within the 
P^C cell the microscopic fluxes are coarse-grained in 



time and space jl9| to supply boundary conditions for 
the continuum domain. 

The spatial decomposition used for the present set-up 
is shown in Fig. ^ The molecular dynamics domain 
ranges from the atomistic wall at z ~ and extends 
to 2: = Icp — 19(7. The continuum fluid dynamics do- 
main comprises Ipc < z < L^, where Ipc — 14.5cr is 
the z— coordinate of the P^C interface and = 50(7 is 
the extent of the whole simulation domain. The center 
of the P^C cell is located aX z — 13.4tT; it has a volume 
Vpc = ^zpcA where A ~ y. Ly and Azpc — 2.2tT 
is the extension along the z direction. The C^P cell is 
placed at a distance 2.2ct from the end of the P^C cell 
and covers a region of Azcp ~ 2.2(T, from z ~ 16.8(7 to 
z = Icp — 19(7. 

In what follows we outline the coupling protocol and 
provide the numerical details used in the present imple- 
mentation. The C— >P coupling represents the most com- 
plicated part of the hybrid scheme; a more detailed ex- 
planation of the method in the frame of the general case 
of unsteady flows with mass, momentum and energy ex- 
changes can be found in reference [23|. The steady flow 
considered here only carries momentum along the x direc- 
tion. Although the mean flux of mass and energy across 
the C and P interfaces is zero, fluctuations in the particle 
system produce perturbative mass currents along the z 
direction which need to be taken into account. This part 
of the C^P scheme is presented in Appendix 1X1 

We now focus on how the momentum flux is exchanged 
between the C and P domains, starting with a discussion 
of the C^P coupling. For the pure Couette shear flow 
considered here, the momentum flux due to the C flow 
across any z =constant surface is given by 

n = Pk - ?77i, (3) 

where P = P{p, T) is the hydrostatic pressure, 77 is the 
dynamic viscosity and 7 = dvx/ dz is the shear rate. The 
value of the dynamic viscosity was measured in a previous 
pure MD simulation via the standard non-equilibrium 
procedure jil[ll|2l|; for p = 0.8 and T — 1.0 we ob- 
tained ?7 = 1.75 ± 0.04. The stress induced by the C- 
flow in the P domain is given by the local momentum 
flux at the C^P interface. Hep. In order to introduce 
this stress into the molecular dynamics domain we add 
an external force Fgxt — — Hcp^ to those molecules 
within the C— >P cell. At any instant of time, t, this 
force is equally distributed among the Ncp{t) particles 
inside the C^P cell, so the external force per particle 
is ^e^t/Npc = -HcF A/Ncp. Note that this exter- 
nal force has a component normal to the C— >P interface, 
which provides the hydrostatic pressure, and a tangential 
component providing the shear stress. The molecules are 
free to enter or leave the C— ^P region, so the number of 
molecules within this region, Ncp{t), and the value of 
the overall external force fluctuate in time. The average 
"pressure force" per particle is P{p,T)/ {ANcp), where 
Ncp ~ 2000 is the mean number of particles within the 
C^P cefl (see Appendix|XI), A = x Ly = 1286cr^ and 
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the pressure P{p,T) is given from the equation of state 
provided by Hess et al. P(0.8, 1) ~ 6.5e/cr3. Such 

a force prevents the escape of particles and, although it 
induces some ripples on the density profile over the C— s-P 
cell, it maintains the correct value of the density along 
the inner part of the MD domain, as seen in Fig. and 
discussed further in the Appendix 1X1 

The shear force is distributed over the particles in the 
same way as described above for the pressure force. In 
this case, the flux of x-momentum to be injected in the 
particle system is rj'jcPi where jcp is the local shear rate 
of the C-flow measured at the C— >P interface. 

We next discuss the P^C coupling. The continuum 
domain is a coarse-grained description of the fluid, there- 
fore any information transferred from the molecular to 
the continuum system needs to be averaged in space 
and time. These averages need to be local in the con- 
tinuum space and time coordinates. To that end, the 
particle quantities are averaged within the P^C cell and 
over a time interval Atav It is important to stress that 
within the P— >C cell each particle's dynamics are not di- 
rectly modified by any external artifact, in other words 
the motion of each particle is uniquely determined by 
the usual molecular dynamics scheme. To ensure consis- 
tency within the hybrid scheme, Atav and the volume of 
the P— >C cell are restricted _26j. For the steady flow cm- 
ployed in this study, the most compelling condition is to 
guarantee that the signal-to-noise ratio of the momentum 
flux is larger than one and for that reason Atav needs to 
increase as 7 decreases We used Atav = lOOr for 
7 < 0.001 and reduced it gradually to lOr for the fastest 
flows considered. 

To solve the equations of motion in the continuum do- 
main we used the finite volume formulation because 
it matches by construction the fiuxes across cells. Since 
the solvent fiow is isothermal, incompressible and there 
is a uniform pressure, the mean x-velocity is governed 
by dvx/dt = vd^Vx/dz^, where v — i]/ p is the kinematic 
viscosity and Vx is the velocity in the x direction. At 
the top of the simulation cell we impose a smooth wall 
in the CFD sense. This wall moves at a constant veloc- 
ity Vx{Lz,t) — Uviaii which creates the shear flow in the 
simulation. The protocol for the P— >C coupling estab- 
lishes the boundary condition for the continuum domain 
at the P^C interface, z — Ipc- The coarse-grained mi- 
croscopic flux of a;-momentum across the P— >C interface, 
whose expression is given in Ref. , is set equal to the 
corresponding value for the C flow at the z = Ipc bound- 
ary, 777PC, where 7pc ~ {vx{lpc + Az) - Vx(lpc))/ Az. 
This condition gives the desired velocity to be imposed 
at the boundary Vx{lpc)- The continuity of velocity is 
ensured by adding a relaxing term in the flux equation 
which drives the C-velocity at the interface towards the 
corresponding averaged P-velocity (see references 
for details). 



III. RESULTS AND DISCUSSION 

In this section we compare the conformational behav- 
ior of the polymers from the MD and hybrid simulations; 
we use two independent MD simulations for comparison. 
We study two MD systems because it is well known 
that two simulations, or experiments, on a tethered poly- 
mer may exhibit considerable variation in conformational 
behavior, even at rather high shear rates. We conclude 
this section with a discussion of the computational costs 
and benefits of the hybrid and classical MD techniques. 

In Figure 121 we show the mean-square end-to-end dis- 
tance E? of the polymer in each of the x, y and z direc- 
tions. This is a standard measure of polymer conforma- 
tion , and its values are related to the values of the 
radius of gyration. At low shear rates the polymer con- 
formation calculated from the hybrid simulation is well 
within the measured conformations of the MD simula- 
tions. At the highest shear rates. Fig. |2Ia) shows the 
conformation of the hybrid polymer to be about 10% 
larger than the polymer in the MD simulation; this dif- 
ference is within the standard deviation of -R^^., which 
is about 15%. Figs|2Ib) and (c) show the y and z com- 
ponents of as a function of shear rate 7. In both 
cases the two MD simulations serve as good indicators of 
the variability of the conformational behavior of a single 
polymer; the conformation of the polymer in the hybrid 
simulation is well within the variations found in the two 
MD simulations, at all shear rates. 

Figure 13 shows the probability of the maximum ex- 
tension from as a function of distance along x, y and z 
directions, for a shear rate of 7 = O.OOlr"^. It is clear 
that the variation of the distributions obtained with the 
hybrid simulation is well within the distribution of the 
two MD simulations. This indicates that not only is the 
average conformation comparable between the two simu- 
lation techniques, but that the probability distributions 
also compare favorably. 

In Figure ^ we show the density, of both the solution 
monomers and polymer, p = N/V as a function of dis- 
tance from the wall for both the MD and hybrid simula- 
tions. The density is calculated in slices of approximately 
O.Olfj perpendicular to the wall. The regular spacing of 
the wall monomers, as two monolayers of a (1, 1, 1) face 
of an fee crystal, induces an ordering in the fluid; this 
ordering is well established [l5l | and persists for approx- 
imately 5tT. At the wall the monomer density variations 
are identical for both the MD and hybrid simulations as 
seen in Fig^Ja). In Fig. EJb) we see that the density of 
both simulations remains the same until the monomers in 
the hybrid system feel the effects of the constant pressure 
condition imposed on the overlap region. The constant 
pressure is implemented as a simple normal force per par- 
ticle on all monomers in the C— >P regime, as discussed 
in the previous section. This force induces a local order- 
ing in the monomers, which in turn creates density fluc- 
tuations. It is noteworthy, however, that these density 
oscillations are much lower than at the atomistic wall. 
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shown for comparison. More recent work on the hybrid 
scheme has estabhshed that we can reduce these density 
fluctations even further, as discussed in the Appendix 
IXI In the MD simulation, the upper waU is identical to 
the lower one, and thus the density fluctuations near the 
former are the same as those at the latter. 

Fig. shows the probability of finding any polymer 
bead in a plane, where the plane slices are 0.2a in thick- 
ness. The two-dimensional probabilities were calculated 
in an analogous way to the one-dimensional probabilities 
discussed above. The shear rate shown is 7 = O.OOI. In- 
spection of the two-dimensional bead distributions indi- 
cates that below a distance of ^ 5a to the wall, the beads 
tend to be ordered in layers parallel to the wall plane. 
This result is not only a consequence of the polymer-wall 
interaction but also an effect of the interaction with the 
solvent. Near the wall the solvent is ordered in layers, as 
in Fig. EJa), and the polymer minimizes the monomer- 
solvent potential energy by adapting its distribution to 
match the locations of the solvent layers. The order in- 
duced by the wall in the polymer structure can be noticed 
even in the isovalues of the probability distribution along 
the wall plane x — y, shown in Fig. I^b), and along the 
z — y plane in Fig. EJc). Over a distance of ~ 6a around 
the attachment position the isovalues of the probability 
distribution in the x — y plane delineate the minimum en- 
ergy lines of the wall atoms LJ potential. In this model, 
the size of the wall atoms was chosen to be the same as 
those of the monomers and solvent particles la. In view 
of Fig. Elb), one should expect that the structure of the 
polymer is quite sensitive to any modification in the de- 
tails of the wall-fluid interaction, owing to cither changes 
in the size of the wall atoms or details of the interaction 
potential. 

In Fig. El we present a comparison of the radius of 
gyration Rg, as calculated from MD simulations in this 
work and that of Aust, Kroger and Hess (ij (AKR) 
who studied a single free polymer in a bath of solvent 
molecules, at a variety of imposed shear rates. The po- 
tential used to describe the polymer and solvent were the 
same in both AKR's work and ours; however AKR used 
a slightly higher density, p = 0.85cr~^, compared to our 
value oi p — 0.8a~^. The simulation of AKR used no 
walls, so the polymer was free to respond to the imposed 
shear so as to best lower the free energy of the system. 
Hence the usefulness of the comparison lies primarily in 
exploring the effect of the wall on the polymer. We see 
that at extremely low shear rates the values of the radii 
of gyration are quite comparable. As the shear rate in- 
creases the value of Rg that we calculate becomes much 
larger than for the equivalent free polymer. This is due 
entirely to the attraction the polymer has with the wall. 

Our results are in agreement with the experimental 
findings of Doyle et al. for individual tethered DNA 
chains under shear flow. For a quantitative comparison 
with these experimental data we evaluated the Weissem- 
berg number, Wi, defined as the product of the shear 
rate and the longest relaxation time of the polymer, i.e. 



the relaxation time at zero shear rate tq. We calculate 
To from the autocorrelation of the polymer extension at 
7 = and obtain tq ~ 2000r. Also, the fractional exten- 
sion is calculated by normalizing the polymer extension 
with its contour length: 0.965 x (TV — 1), where iV = 60 
is the the number of monomers and 0.965(j is the mean 
separation between two consecutive beads Using 
this value of tq we plot in Fig. the normalized mean 
fractional extension along the flow direction versus the 
Weissemberg number, along with the the expermimental 
results of Doyle et al. The results obtained with both, 
the MD and hybrid simulations are in very good agree- 
ment with the experimental data for the range of shear 
rates considered here. 



Figure |H1 shows the end-to-end volume of the polymer, 
measured as the product of the three components of the 
end-to-end vector. This quantity gives an estimate of the 
space that the polymer explores during its motion. This 
volume increases for increasing shear rate and reaches a 
maximum value around 7 ^ 0.002t~^. As the shear rate 
is further increased the volume accessible to the polymer 
decreases monotonically. This behavior of the end-to-end 
volume is quite similar to the findings of Doyle et al. 
concerning the amplitude of the fluctuation of the chain 
extension. As the shear rate was varied, they found that 
fluctuations reached a maximum size at Wi ~ 5.1. Using 
the estimate tq ~ 2000t, we find that the maximum end- 
to-end volume occurs at about Wi ^ 4. In fact, the size 
of the fiuctuations is determined by the magnitude of 
the volume made available by the polymer motion; or, 
in other words, larger fluctuations increase the explored 
volume. We shall present a more detailed comparison 
with the results of Doyle et al. [6| in future work. 



There are computational costs to the hybrid method 
that are not present for classical MD. These include simu- 
lating the continuum regime and the calculations arising 
from the coupling procedure within the overlap region, 
e.g. particle insertion and deletion and the evaluation of 
the particle stress tensor. For the flow considered here, 
the solution of the continuum flow required around 0.01% 
the time needed for a LJ force calculation. In general, the 
computational time spent in simulating the continuum 
region depends on the problem considered, but in any 
case it will always be much smaller than the MD force 
evaluation for the solvent. Furthermore, the calculation 
of the C-flow occurs once for every ~ 20 LJ force cal- 
culations, which ensures extra savings in computational 
time. As shown in Appendix ^ the coupling protocol, 
within the overlap region, is very efficient: only 0.01% of 
the total computational time was spent in particle inser- 
tion and deletion while around 5% in the evaluation of 
the particle stress tensor. The hybrid code as tested here 
needs less than half the solvent particles, thus the overall 
savings in computational time is considerable. 
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IV. CONCLUSION 

In this paper we have compared a newly developed 
hybrid MD/CFD code to a traditional MD simulation for 
a single polymer tethered to a wall undergoing shear flow 
in Couette geometry. We find that the two methods give 
comparable results for the conformation of the polymer 
within measured uncertainty. 

Our results indicate that the coupling protocol of the 
hybrid code requires around 5% of the computer time 
compared to the Lennard- Jones part of the code. Most 
of the CPU time devoted in the "coupling" protocol is 
spent in the evaluation of the particle momentum flux; 
while insertion and extraction of particles is rather fast, 
taking less than 1% of the overall CPU time. 

This implies that, compared with a traditional MD 
simulation, the amount of computational time saved by 
the hybrid scheme is proportional to the volume of the 
simulation that is described by the coarse-grained model 
(CFD). In traditional MD simulations of interfacial phe- 
nomena finite size effects significantly alter the local in- 
terfacial dynamics, and they can only be reduced by in- 
creasing the volume of the simulation box that surrounds 
the interfacial region of interest. This means that most 
of the computational cost is likely to be spent in the res- 
olution of the bulk flow. In this paper we have shown 
that this drawback disappears when using a proper hy- 
brid MD-CFD scheme. To that end, we considered a 
problem which is very sensitive to small changes in the 
surrounding fluid environment: the motion of a single 
tethered polymer under shear flow. The excellent agree- 
ment found in the comparisons with the full MD results 
indicates that the hybrid scheme indeed eliminates finite 
size effects even in relatively small systems. This means 
that hybrid simulations can be expected to significantly 
reduce the computational cost of appropriate interfacial 
problems. 

Apart from the savings in CPU time, the hybrid 
scheme enables to us gather information from all rele- 
vant time and length scales, so it is well suited to treat 
multiscale problems where bulk fluid flow plays an im- 
portant role; other examples include crystal growth from 
fluid phases, wetting phenomena and membrane dynam- 
ics under flow. 

We regard the results of the present work as an en- 
couraging sign for future simulations, and we plan to ex- 
plore various selected interfacial systems in forthcoming 
research. 
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APPENDIX A: MASS, LONGITUDINAL 
MOMENTUM AND ENERGY FLUCTUATIONS 

In this work the mean solvent flow carries no longitudi- 
nal momentum along the z-direction and has a constant 
mean density. However, we observe that the polymer 
motion induces density and longitudinal velocity fluctu- 
ations within the particle region that induce currents of 
mass and longitudinal momentum travelling along the 
simulation box. These perturbative currents have to be 
controlled at the C^P interface. We need to ensure 
that the mean mass flux across the z — Icp interface 
is zero, but in such a way that any pressure waves leave 
the simulation box once they reach the C^P interface. 
In other words, we need to prevent any pressure waves 
from bouncing back at the C^P interface in the MD 
region. 

The average number of particles crossing the C^P in- 
terface per unit time is given by Ncp = A{pvz)cp- Zero 
mass flux is ensured by equating this rate Npc to the 
rate of insertion of molecules into the particle system 

^1 . In the calculations presented here we used an- 
other control equation which provides a finer control on 
the particle density near the C^P interface. This ap- 
proach is based on relaxing the local density at the C^P 
buffer to a prespecified value po, 

Ncp^ — {{p)cp-po) (Al) 

where Vcp is the volume of the C^P cell, {p)cp is its 
local the particle density averaged over ISiav and is a 
relaxation time which controls the rate at which the den- 
sity fluctuations within C^P cell are smoothed out. We 
set the value of r„i slightly smaller than the time needed 
by a sound wave to cross the C^P cell (~ 0(1)t). This 
procedure ensures that fluctuations carrying mass and 
longitudinal currents are damped at the C^P cell and 
do not bounce back to the inner part of the MD do- 
main. According to Eq. ljAl|l . particles are extracted 
if Npc < and, as explained in j[2Q], the first particles 
to be extracted are those closer to the C^P interface. 
If Npc > 0, new particles are inserted with a velocity 
extracted from a Maxwellian distribution with mean ve- 
locity Vy — Vz = and Vx — jz and temperature T = 1.0. 
The insertion of particles in liquids is not a trivial task, 
however, and it is addressed by the usher algorithm for 
particle insertion The value of po in Ea. (|Al|l was 

set to a slightly smaller value, po = 0.65, than the mean 
density 0.8. This reason for this choice is, first, to allevi- 
ate the computational cost of insertion (see Appendix IB|) 
and, second to reduce the amplitude of the ripples of the 
density profile at the C^P buffer, as shown below. For 
a liquid with p ~ 0.8 the usher algorithm needs around 
30 iterations to insert a LJ atom at a location where the 
potential energy equals the mean specific potential en- 
ergy of the system |l9j , where each iteration corresponds 
to the evaluation of a single-particle force. If the density 
is decreased to 0.65, it only needs about 15 iterations. 
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Figure 6 compares the density profile resulting from us- 
ing p — 0.65 in Eq. I|A1() with that arising from a pure 
MD simulation. The "hybrid" density profile presents 
some ripples whose amplitude is damped after around 
3(7, whereas inside the P^C cell the hybrid density pro- 
file perfectly matches the density within the bulk. 

As long as the fiuid is isothermal and there are no 
mean pressure gradients, the mean energy fiux across the 
C-^P interface is zero. We therefore only need to guar- 
antee that the specific energy of the newly inserted par- 
ticles matches that of the ensemble. The kinetic energy 
is matched by inserting new particles with a Maxwellian 
distribution, as stated above. In order to match the po- 
tential energy, new particles are inserted at sites where 
the inter-particle potential energy equals the chemical 
potential of the system, thereby ensuring the Widom in- 
sertion criterion. 



APPENDIX B: COMPUTATIONAL COST OF 
THE HYBRID METHOD 

We compare the computational cost of the coupling 
subroutines with those pertaining to the MD part of the 
hybrid scheme. This comparison was made using the 
gprof command available in the package of the f 77 com- 
piler. One of the parts of the hybrid scheme for which one 
may expect a certain cost in computational time is parti- 
cle insertion. Table 1 presents some results obtained for 
different shear rates and values of the density po in Eq. 
(|IH) . Typically, Eq. requires around 5 insertions 

per time interval r and around 15 iterations per parti- 
cle (each interation involving a single- force evaluation). 
Therefore, for a time step of Aip = 0.0075r, the inser- 
tion of new particles needs typically about one extra force 
evaluation per time step. This number is very small when 
compared with the number of force evaluations needed in 
the MD system, which is on the order of the number of 
particles Np ~ 10**. This estimate is consistent with our 
findings concerning the computational cost. As shown in 
table 1, in hybrid calculations using po = 0.65, the time 



spent in the insertion/extraction subroutines was about 
1.5 X 10"'' times the time spent in the force evaluation 
and around 0.9x 10^^ if one includes the Verlet list eval- 
uation. This performance confirms the extremely high 
efficiency of the usher algorithm for particle insertion. 

As a matter of fact, the dominant cost of the hybrid 
scheme resides in the evaluation of the particle momen- 
tum flux. Its cost in CPU time was about 0.06 times the 
cost of the force plus Verlet list subroutines. We note 
that the implementation of this part of our code was not 
constructed in an efficient way because we evaluated the 
particle momentum flux at each MD time step. Consider- 
ing that for the evaluation of (jp) we used measurements 
of jp separated by its decorrelation time, about 0.06t |^, 
we could have measured the particle flux roughly every 10 
time steps and further reduced that ratio by a factor 10. 



7(r-i) 


po 








CPU[insErt\ 
CPU\force] 


0.001 


0.8 


3.68 


25.4 


0.015 


2.710"* 


0.010 


0.8 


3.64 


25.9 


0.015 




0.010 


0.65 


8.25 


14.7 


0.006 


0.910"* 


0.005 


0.65 


4.34 


16.3 


0.010 





TABLE I: Details of the particle insertion in several hybrid 
simulations done at shear rate 7 . Using po in Eq. ijAl^ . 
the average rate of particle insertion was Nin and the aver- 
age number of iterations needed by the usher algorithm to 
insert a new particle Ee is the relative error in 

the energy upon insertion (the relative difference between the 
target potential energy and the potential energy at the inser- 
tion site). In the last column we show the ratio between the 
CPU time used by the insertion/extraction subroutines and 
the CPU time used by the force subroutine plus the Verlet 
neighbor list. 



in the continuum domain was very small compared with 
the MD force subroutine, by a factor of less than 10"*. 
In general, the computational time required to solve the 
continuum system will, of course, depend on the specific 
problem solved. 
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FIG. 1: The domain decomposition of the hybrid scheme. 
The polymer is embedded within the particle region (P) which 
is described by molecular dynamics, including the atomistic 
(lower) wall and the solvent (Lennard- Jones particles). Fluid 
flow within the continuum region (C) is described by an un- 
steady Stokes equation and is solved using finite volumes. The 
handshaking region contains the C^P and the P— >C cell, 
where the two-way exchange of information is established. 
The P and C domains overlap within Ipc < z < Icp- The 
area of the P^C and C— »P cells is the surface of the sys- 
tem in the periodic directions, A — L^Ly. The Couette flow 
moves along the x direction driven by the velocity imposed 
by the upper boundary condition, which corresponds to the 
upper wall velocity in the pure MD simulations, u^aii- The 
magnitudes of each length shown in the figure are given in 
Sec. EO 
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FIG. 2: The x, y and z components of the mean square end- 
to-end vector, R^, are shown as a function of shear rate for 
two independent MD simulations and one hybrid. Error bars, 
not shown, are approximately 15%. The cc-component of 
increases as the shear rate increases, while the y and z com- 
ponents decrease. At low shear rate the hybrid simulation is 
well within the variation of the MD simulations. At the high- 
est shear rate the values for R^^ agree within the measured 
uncertainty. 




FIG. 3: Probability of finding a monomer in the x (a), y (b) 
and z (c) coordinates in a flow with shear rate 7 = O.OOlr"^. 
Comparison is made between the result obtained with the 

hybrid scheme and the outcome of two pure MD simulations 
with different initial conditions. 
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FIG. 4: Monomer and solvent density as a function of distance 
(z) from the wall, for shear rate 7 = 0.005t~^. (a) shows 
the density fluctuations near the lower wall. In (b) dashed 
lines indicate the locations of the coupling buffers used in the 
hybrid scheme, P^C and C^P. Outside the C — >P region, 
the monomer density is unaffected by the hybrid scheme 
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FIG . 5 : Probability of finding a monomer inthe x — z (a), x — y 
(b) and y — z (c) planes in a fiow with shear rate 7 = O.OOlr"^ . 
The maximum of the probability distribution is located near 
the attachment site. The shaded region corresponds to an 
iso-probability value of 0.021 and the values of consecutive 
iso-probability contour lines are separated by 0.01. The his- 
tograms were obtained from the calculation of a pure MD 
simulation with a total simulation time of 78750r. 




FIG. 6: Comparison of the radius of gyration of the polymer 
as calculated by MD simulation in this paper, and that of 
Aust et al. [Tlj . The higher Rg found in this work is due to 
the attraction between the wall and the polymer. 
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FIG. 7: The fractional elongation along the flow direction 
versus the Weissemberg number, Wi=ro7. The longest decay 
time at zero shear rate obtained from our data is tq — 2000r. 
Comparison is made with the experimental results of Doyle 
et al. 01 on tethered polymers. 
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FIG. 8: The end-to-end volume of the polymer as a function 
of the shear rate. The end-to-end volume is defined by the 
product of the components of the end-to-end vector, {R^ x 
Rl X Rlfl^. 



